clear
version 14.0
set more off

/*
* Translate shapefile to stata format
cd "../maps/"
shp2dta using my_county_map_48states_albers, ///
		database(county_db_usa) coordinates(county_coord_usa) genid(id) replace

shp2dta using my_state_map_48states_albers, ///
		database(state_db_usa) coordinates(state_coord_usa) genid(id) replace		
		
shp2dta using county_my_region2, ///
		database(county_db_my_region2) coordinates(county_coord_my_region2) genid(id) replace		

shp2dta using state_my_region2, ///
		database(state_db_my_region2) coordinates(state_coord_my_region2) genid(id) replace			
	
		*/



* Specify the year of cash rent data to use in the regression
local reg_year 2013

use "..\dataAnalysis\rent_climate_data", clear



rename stfips stfip
rename cofips cofip
merge 1:1 stfip cofip using "../maps/county_db_my_region2", keep(match using) nogen

* Drop counties in SW New Mexixo and Texas to allow room for the legend
drop if stcofips==35023 | stcofips==35017 | stcofips==35029 | stcofips==35013 | stcofips==48141 | stcofips==35035 ///
		 | stcofips==35051 | stcofips==48229 
		 
gen my_region=(lrrsym=="M" | lrrsym=="H" | lrrsym=="F" | lrrsym=="G" | lrrsym=="O" | lrrsym=="L" ) 
* Drop New York
replace my_region=0 if stfips==36		 

summ cash_rent_non2013, detail
format cash_rent_non2013 %9.0f
gen map_rent=.
replace map_rent=7 if cash_rent_non2013>250 & cash_rent_non2013!=.
replace map_rent=6 if cash_rent_non2013>200 & cash_rent_non2013<=250
replace map_rent=5 if cash_rent_non2013>150 & cash_rent_non2013<=200
replace map_rent=4 if cash_rent_non2013>100 & cash_rent_non2013<=150
replace map_rent=3 if cash_rent_non2013>50 & cash_rent_non2013<=100
replace map_rent=2 if cash_rent_non2013>25 & cash_rent_non2013<=50
replace map_rent=1 if cash_rent_non2013>10 & cash_rent_non2013<=25
replace map_rent=-1 if cash_rent_non2013==. & my_region==1
replace map_rent=-2 if cash_rent_non2013==. & my_region==0

label define rentlab 7 "(250,385]" 6 "(200,250]" 5 "(150,200]" 4 "(100,150]" 3 "(50,100]" 2 "(25,50]" ///
			1 "(10,25]" 0 "Dropped" -1 "" -2 ""
label values map_rent rentlab
tab map_rent

spmap map_rent using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			ndfcolor(gs11) clmethod(unique) legorder(hilo) ///
			legend(position(8) size(3.4))    ///
			fcolor(white white "237 248 233" ///
					"199 233 192" "161 217 155" "116 196 118" ///
					"65 171 93" "35 139 69" "0 90 50" )
/*		
spmap map_rent using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(white) osize(medthick)) ///
			ndfcolor(gs11) clmethod(unique) legorder(hilo) ///
			legend(position(8) size(3.4))    ///
			fcolor(gs11 gs14 "237 248 233" ///
					"199 233 192" "161 217 155" "116 196 118" ///
					"65 171 93" "35 139 69" "0 90 50" )
*/
graph export "../figures/rent_map.png", replace	width(1200)


summ water_deficit, detail
format water_deficit %9.0f
spmap water_deficit using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(YlOrRd) ndfcolor(white) ndlabel(" ")   ///
			clm(custom)  clbreaks(2 5 7 10 15 20 25 35) ///
			legend(position(8) size(3.5)) title("Water Deficit") ///
			plotregion(margin(l=7 r=2 t-5)) name(water_deficit, replace)
graph export "../figures/water_deficit_map.png", replace	width(1200)	






summ water_surplus, detail
format water_surplus %9.1f
spmap water_surplus using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(GnBu) ndfcolor(white) ndlabel(" ")   ///
			clm(custom)  clbreaks(0 0.1 0.3 0.5 1 1.5 3.0 5.0) ///
			legend(position(8) size(3.5)) title("Water Surplus") ///
			plotregion(margin(l=7 r=2 t-5))  name(water_surplus, replace)
graph export "../figures/water_surplus_map.png", replace	width(1200)				



summ gdd, detail
format gdd %9.0f
spmap gdd using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Oranges) ndfcolor(white) ndlabel(" ")   ///
			clm(custom)  clbreaks(9 12 15 18 21 24 27) ///
			legend(position(8) size(3.5)) title("GDDs") ///
			plotregion(margin(l=7 r=2 t-5))  name(gdd, replace)
graph export "../figures/gdd_map.png", replace	width(1200)
		



summ dday34C, detail
format dday34C %9.0f	
spmap dday34C using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Reds) ndfcolor(white) ndlabel(" ")  ///
			clm(custom)  clbreaks(0 1 2 4 10 20 40 61) ///
			legend(position(8) size(3.5)) title("EDDs") ///
			plotregion(margin(l=7 r=2 t-5)) name(dday34C, replace)
graph export "../figures/dday34C_map.png", replace	width(1200)


*------------------------------------
* Combine graphs
*-------------------------------------
graph combine water_deficit water_surplus gdd dday34C, ///
		cols(2) imargin(vsmall) iscale(*1.27) graphregion(color(white)) altshrink ysize(4) xsize(5)
*graph combine perc_chg_rent_total perc_chg_rent_deficit perc_chg_rent_surplus perc_chg_rent_heat, ///
*		cols(4) imargin(vsmall) iscale(*1.27) graphregion(color(white)) altshrink ysize(2.25) xsize(7.5)

graph export "..\figures\hist_climate_combined.png", replace	width(1600)	


*------------------------------------
* Irrigated cash rent
*-------------------------------------
summ cash_rent_irr2013, detail
format cash_rent_irr2013 %9.0f
spmap cash_rent_irr2013 using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Greens) ndfcolor(white) ndlabel(" ")   ///
			clm(quantile) cln(7) ///
			legend(position(8) size(3.4)) 
graph export "../figures/rent_irr_map.png", replace	width(1200)	
	
*------------------------------------
* Irrigation, Precip, and water holding capacity
*-------------------------------------
gen proportion_irrigated = acresh_Irr2012 / (acresh_Irr2012 + croplandNon_acres)
summ proportion_irrigated, detail
format proportion_irrigated %9.2f
spmap proportion_irrigated using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(YlGnBu) ndfcolor(white) ndlabel(" ")   ///
			clm(custom)  clbreaks(0 0.05 .15 .25 .35 .45 .57) /// 
			legend(position(8) size(3.4)) 
graph export "../figures/irrigated_map.png", replace	width(1200)	

summ ppt, detail
format ppt %9.0f
spmap ppt using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Blues) ndfcolor(white) ndlabel(" ")   ///
			clm(quantile) cln(7) ///
			legend(position(8) size(3.4)) 
graph export "../figures/ppt_map.png", replace	width(1200)	


summ rootznaws, detail
format rootznaws %9.0f
spmap rootznaws using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Purples) ndfcolor(white) ndlabel(" ")   ///
			clm(quantile) cln(7) ///
			legend(position(8) size(3.4)) 
graph export "../figures/rootznaws_map.png", replace	width(1200)	


summ winter_wheat, detail
format winter_wheat %9.2f
replace winter_wheat=0.99 if winter_wheat>1 & winter_wheat!=.
spmap winter_wheat using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Blues2) ndfcolor(white) ndlabel(" ")   ///
			clm(custom)  clbreaks(0 0.1 0.2 0.3 0.4 0.5 0.6 1) /// 
			legend(position(8) size(3.4)) 
graph export "../figures/winter_wheat_map.png", replace	width(1200)	
			
capture drop my_wheat			
gen my_wheat=0
replace my_wheat=1 if winter_wheat>0.1
replace my_wheat=. if winter_wheat==.
spmap my_wheat using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Set2) ndfcolor(white) ndlabel(" ")   ///
			clm(unique)   /// 
			legend(position(8) size(3.4)) 			
			
			
summ text_sandy text_mod_sandy text_medium text_mod_clay text_clay
gen my_sandy=text_sandy + text_mod_sandy
gen my_clay=text_mod_clay + text_clay
gen my_silt=text_medium
summ my_sandy my_clay my_silt

spmap my_sandy using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Purples) ndfcolor(white) ndlabel(" ")   ///
			clm(quantile) cln(7) ///
			legend(position(8) size(3.4)) 
			
spmap my_clay using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Purples) ndfcolor(white) ndlabel(" ")   ///
			clm(quantile) cln(7) ///
			legend(position(8) size(3.4)) 
			
spmap my_silt using "../maps/county_coord_my_region2", id(id) ocolor(none ..) ndocolor(none ..)  	///
			polygon(data("../maps/state_coord_my_region2") ocolor(gs14) osize(medthick)) ///
			fcolor(Purples) ndfcolor(white) ndlabel(" ")   ///
			clm(quantile) cln(7) ///
			legend(position(8) size(3.4)) 
